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An atomic-scale model for silicate solutions is introduced for investigation of the nucleation process 
during zeolite synthesis in the absence of a structure directing agent. Monte Carlo schemes are 
developed to determine the equilibrium distribution of silicate cluster sizes within the context of this 
model. How the nucleation barrier and critical cluster size change with Si monomer concentration 
is discussed. Distance and angle histograms as well as ring size distributions are calculated and 
compared with known zeolite structures. The free energies of critical clusters are compared with 
those for small clusters of a-quartz. 
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I. INTRODUCTION 



Zeolites are CEMStalline microporous materials that have found a wide range of use as catalysts, molecular sieves, and 
ion exchangers. Elu A typical example is ZSM-5, used as a cracking catalyst in the refinement of crude oil. Essentially 
all refinement of petroleum occurs through fluid catalytic cracking with a H-ZSM-5 zeolite catalyst. Essentially all 
reforming occurs with zeolitic catalysts, typically H-ZSM-5 or faujasite-type structures. A significant fraction of air 
separation into N2 and O2 occurs via a process involving the zeolite Li-X. Roughly 25% by weight of the powdered 
laundry used in washing machines is the zeolite Na-A, which softens the water for low-phosphate detergents. Finally, 
zeolites are widely used in remediation applications. For example, victims of the three-mile-island accident were fed 
ion-exchange zeolites to remove radioactive traces of ^^^Cs+ and ^"Sr^^. 

Zeolites continue to be synthesized at a furious pace. The major current limitation of our ability to synthesize .geolites 
with tailored properties is an incomplete understandiag of the fundamental steps that occur in zeolite synthesisjS which 
is normally carried out at hydrothermal conditions.u Important factors that determine js¥rtich zeolite will be made 
include the Si/Al ratio, the alkalinity, the cation species, and the presence of templates-clcl A significant amount of 
experimentation, including NMR spectroscyaay.-X.xay diffraction, and neutron scattering, has been done in an effort to 
characterize the growth kinetics of zeolites ,aQi3'ElBE2lO't3 yet the mechanism of nucleation is still not fully understood. 
One of the difficulties arises from the accessible length scale in experiments. NMR is able to provide information on 
the 0.5-1.0 nm scale. Diffraction, on the other hand, is limited to providing detailed information on the 5 nm and 
above length scale. The nucleation process falls, unfortunately, within this gap and is difficult to beidentified directly. 
High temperature calorimetry provides thermodynamic data but no direct structural information.EJ Recently, it has 
been possible to access the nucleation scale with a combination of electron microscopy. X-ray, and neutron diffraction 
methods. □ A working hypothesis formed is that, at least for template-mediated synthesis of ZSM-5, the nucleation 
event involves roughly 8 disordered unit cell precursors aggregating into a "primary unit" . This aggregate forms the 
nucleation core, to which additional units are added. It is unclear whether the larger 5-10 nm globule is already 
crystalline or is amorphous. As the globule grows, a small crystallite begins to form. This type of nucleation on the 
length scale of a few unit cells is consistent with typical nucleation behavior of other electrolytes in aqueous media. 
In particular, the primary unit size for ZSM-5 is roughly 2 unit cells on a side, consistent with the nucleation scale 
of 2-3 unit cells for typical electrolytes condensing in aqueous media. Further support is provided by the observed 
primary unit size of 2.6 gm for zeolite beta (unit cell of 1.3 x 1.3 x 2.4 nm)L3 and 1.6 nm for zeolite SOD (unit cell 
of 0.9 X 0.9 X 0.9 nm).u These primary units were not detected under conditions for which no zeolite was formed 
(i.e., without template). The concentration of primary units decreased as the initial stages of nucleation and growth 
took place. All of these data show the strong correlation between the presence of the primary units and the rate of 
nucleation. This hypothetical mechanism for ZSM-5 synthesis, however, has not been unambiguously demonstrated 
to be correct. 

The hypothetical zeolite nucleation event involves few enough atoms that it may be investigated computationally. 
Moleculaii.di"namics has been used to examine the early stage of the polymerization process of silicate ions in aqueous 
solution.Ejllj This dynamical approach yields important information of the reaction mechanism, but is limited to the 
case where the system is well above the saturation concentration and the free eneiigjtJparrier is low. Methods that 
can investigate the nucleation event for simple solids have prpiitly been developed .113113 Related methods have been 
developed in the context of ion-induced aerosol nucleationJljcj These methods make use of transition state theory, 
and so they avoid the time-scale problems associated with a direct simulation of the nucleation event. 

As a first step towards fundamental understanding of the zeolite nucleation process, we present here a Monte Carlo 
study of silicate solutions, using a model with explicit Si and O atoms. We aim at constructing an equilibrium 
distribution of clusters that provides the free energy barrier of nucleation. Novel reactive moves that alter the 
connectivities of the silicate clusters, as well as other Monte Carlo moves, are used to equilibrate the system. No a 
priori building blocks are assumed in the simulation. Experimentally known features of aqueous silicate solutions are 
used to calibrate the model. We analyze the structure of the clusters found in the silicate nucleation process, and we 
compare the results with known zeolite structures. 

This article is organized as follows. We describe the model for the silicate solution in Sec. II. We describe in Sec. fll 
how to compute the equilibrium distribution of clusters under different experimental conditions. The Monte Carlo 
simulation schemes are described in Sec. IV. Results are given in Sec. and discussions are made in Sec. Vl. Finally, 



we conclude in Sec. VII 



II. MODELING THE SILICATE SOLUTION 



A key component of simulation of silicate nucleation is a proper forcefield to describe the energetics of the silicate 
species, hydrogel, cations, and templates. We focus on the pure-silica case for the nucleation process, and model the 
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silicate clusters by explicit, charged Si and O atoms. It is advantageous to distinguish the atoms by their connectivities. 
Each Si atom is sp^ bonded to four O atoms. An O atom that is bonded to one O atom is denoted by O*^^-', and an O 
atom that is bonded to two Si atoms is denoted by O*^^-'. We denote by Si'-*^ a Si atom that is bonded to i O^^-* atoms 
and (4-j) O*-^-* atoms and call this an i-connected Si. Neutralizing charges in the solution are treated implicitly by 
adding a uniform, neutralizing background charge distribution to the system. Hydrogen atoms terminating a SiO*-^-* 
bond are also treated in an implicit fashion. 

Since the nucleation and crystallization of the silicate species involves the formation of covalent bonds, we will 
need to use a-j|eajetiM2|jMoe£eld. Such forcefields have been developed for silicate species, and they continue to 
be developed. E^OS'EjoESEa Here we use a modified version of the potential proposed by Vashishta et al.E3 This 
potential includes pairwise and three-body interactions between explicit atoms. The two-body part accounts for the 
steric repulsion. Coulomb interactions, and charge-dipole interactions. The explicit form of the two-body interaction 
is 

V2 = — ^ + — — ^ exp -r/r4s). (1 

Here Hy and 'qij represent the strengths and exponents of the steric repulsion. The qi and Ui are the effective charge 
and electronic polarizability of the zth ion, respectively. The decay length r4s is taken to be 4.43 A.E3 The expression 
for the three-body interaction is 
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Bjik exp 
0, Ti, > ro or r^k > ro 



{cos9jik - cos 9 jik), nj < ro and r^fc < rg 



(2) 



where Bjik is the strength of the three-body interaction, I is a constant, rg is a cutoff distance, 9jik is the angle 
subtended by the vectors r^j and Vik with the vertex at i, and 9jik is a constant angle. The constants are summarized 
in Table |. 

One of our modifications of the potential takes into account the solvent effects, which are modeled by applying 
a distance-dependent dielectric constant e = 4r in eq. |l|.EII Our second modification attempts to incorporate the 
deprotonation effects. Since g(0(2)) = -0.80 e, as seen in Table ||, each Si-0 bond contributes -0.40 e to the O'-^-' 
charge. Silicate solutions at zeolita-synthesis condition have pH = 11 — 12, at which about half of the Si-O'""^^ • • -H 
bonds of monomers are dissociated.c3 Assuming an O- • -H pair has a total charge of -0.40 e, and the de-hydrogenated 
O^^^ has a total charge of -1.40 e, we assign the average value of -0.90 e to the O^^^ atoms, thus taking into account 
the terminal hydrogens in an implicit way. 



III. DETERMINING THE CLUSTER-SIZE DISTRIBUTION 



Direct simulation of the aqueous silicate solution will fail to observe the nucleation event, because of the extremely 
low probability for the relevant fluctuations to occur. Despite this difficulty, it has been shown that the ppabability 
distribution of clusters can be obtained from simulation of a single cluster in the grand canonical ensemble. EElo Such 
a distribution can be used to quantify the free energy of clusters and to predict the nucleation rate from theorytZlu 

Consider a system specified by volume temperature T, and the chemical potentials fisi for Si and fio for O, 
where fisi and fio are implicitly related to the concentration and pH of the system. The reference system is taken 
to be the monomeric Si(OH)2 02~ ion system. A uniforai, neutralizing background charge density is added to the 
reference system so as to make the system energy finite.tll When the system becomes super-saturated but confined 
in the meta-stable region, polymerization of the silicates leads to larger clusters. We use the following geometrical 
criterion for a cluster:E3 a Si-0 bond is formed when their distance is within 2.60 A, each Si atom is bonded to exactly 
four O atoms, an O atom can be bonded to either one or two Si atoms, no Si-Si bonds or 0-0 bonds can be formed, 
and a cluster is a set of atoms that are connected by bonds. We classify a cluster by its numbers of Si and O atoms, 
denoted by ng; and no, respectively. Since nucleation is a rare event, the population of clusters is very low, jHiid the 
interaction between clusters is negligible. In this case, the partition function of the system can be written aso 

lnS(^Si,MO, V^,T) = ^ exp [/3(nsiMsi + "-0M0)] Qnsi,"o ) (3) 

nsi,"o 

where /3 = l/lksT) is the reciprocal temperature, and the partition function for an nsi, no-cluster, Q„si,no i defined 
as 

V_ 

Asi3"^'A;7^°nsi!7^o!' 
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^„si,no = yrfrsi"^'rfro"°(5(R„n)«^(rsi"^-,ro"°)exp[-/3[/(rsi"^',ro"°)]. (4) 

Here Asi and Aq are the thermal wavelengths for Si and O atoms, respectively, Rem is the center of mass of the 
cluster, and the geometrical cutoff function ^(rsi"^', ro"°) is 1 if rsi"^',ro"° satisfy the cluster criterion, and it 
is otherwise. Note that i(;(rsi"s', ro"") is invariant upon translation of the coordinates. The potential energy 
C^(rsi"^S ro"°) accounts for the intra-cluster interaction. The translational degrees of freedom of the cluster are 
integrated out to get V. Note that the configurational integral in cq. ^ is independent of V as long as V is larger than 
the cluster size. |— ■ 
The average number of clusters of size nsi, no is given byt3 

{N{nsi,no)) = exp [/3(^sinsi + Mo»^o)] Qnsi,no- (5) 
The partition function for a single cluster is given by 

^cluster = ^ exp [/3(nsi/Xsi + noAio)] Qnsi.no- (6) 

That is, the partition function of the system of non-interacting clusters is the exponential of the partition function 
for a single cluster. 

We define Aft as the excess free energy of a cluster over an isolated monomer. It is given by 

/3An{nsuno,^lsi,^J'0,T) = -In ^^^^ . (7) 

A natural choice of the nucleation coordinate is nsi. Therefore, we are mainly concerned with the dependence of the 
free energy on nsi- The free energy of the nsi-cluster, which is the union set of ngj, no-clusters with all possible no, 
is defined as 



/3Ar2(nsi, Aisi, MO, T) = - In < ^ exp [-/3A17(nsi, no, fxsi, Mo, T)] ? 

\ no ) 



(8) 



In the simulation only one cluster is present. The variables that specify the system are ^Si, Mo, and T. As shown in 
Sec. an arbitrary system volume V can be used in simulation, because-the particle insertion/removal moves are 
proposed according to a geometrical criterion instead of the system volume.E3 To simplify the notation in the following 
sections, we define the fugacities zsi — exp(/3/iSi)/Asi^ and zq — exp(/3Mo)/Ao'^- Equation || can be rewritten as 

(jV(nsi,no)) _ ^Si"^-^o"" ^ 

V ~ nsiino! ^ ' 



IV. MONTE CARLO SCHEMES 



We use several different Monte Carlo moves to sample different degrees of freedom in our system. At all times, only 
a single cluster is present in the simulation box. Both the configurational degrees of freedoms of the cluster and the 
numbers of Si and O atoms in the cluster are sampled by the Monte Carlo moves. These moves are chosen randomly 
from a set with predetermined probabilities. The move frequencies are different for small and large clusters so as to 
achieve effective sampling. The cluster criterion is enforced at all times during the simulation to avoid dissociation of 
the cluster. 



A. Reactive Moves 



The covalent bond formed between Si and O atoms has a magnitude that is much larger than thermal fiuctuations 
at the synthesis temperature. Therefore, direct insertion or removal of particles will fail, due to the rough energy land- 
scape. Thus, we developed smart Monte Carlo reactive moves to overcome this problem. The Monte Carlo reactive 
moves are divided into three types, the SiO„ insertion/removal move, where n = 2 or 3, the hydrolysis/condensation 
move, and the prune-and-graft move. We demonstrate in Fig. ^ how the SiO„ insertion/removal move and hydroly- 
sis/condensation move help the system to equilibrate the chemical composition (nsi,no) of the cluster. For example, 
a SiOa insertion move brings the cluster from (nsi = 2, no = 7) to (nsi = 3, no = 10), as shown by the arrow in 
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Fig. The prune-and-graft move does not change nsi and no and thus is not shown on the lattice. We here describe 
the Si02 insertion/removal moves. Other moves are described in Appendix 

A SiO„ insertion move attempts to graft a SiO„ group onto the surface of the cluster. A SiO„ removal move 
attempts to remove a SiO„ group from the surface of the cluster. Figure ^ illustrates the Si02 insertion/removal 
moves. A Si02 insertion move proceeds as follows: 

1. Pick randomly a pair of O^^-* atoms from a Si02 insertion list that is constantly updated during the simulation. 
The two O atoms must satisfy certain geometrical criteria, which eliminate geometrically impossible pairs so 
as to improve efficiency of the moves. Only distinct pairs of atoms are included in this list; that is, only the 
identity of the pair of atoms, and not the order, is included in this list. Consider the two Si atoms to which this 
pair of O*^^-' atoms are bonded. Mark as group 1 the first Si atom and the O*^^^ atoms bonded to it. Mark as 
group 2 the second Si atom and the O*^^^ atoms bonded to it. 

2. Generate trial local moves for group 1 and group 2, according to the O*^^^ connectivity i of each Si atom. When 
i = 1, generate a random orientation for the SiOa group. When ? = 2, rotate the Si02 around the axis passing 
through the two O^^^ atoms bonded to the Si atom. When z = 3, no move is made. In all cases, multiple trial 
moves are generated, each move consisting of a modification of both groups. One of the trial moves is chosen 
with a probability proportional to p(roo) exp(— /3[/)/(47rroo^)- Here roo is the distance between the pair of 
O'^' atoms, and p(roo) droo is the probability that two chosen O atoms in a monomer have a distance between 
roo a-nd roo + droo- This probability is predetermined by a short hybrid Monte Carlo run for a monomer. 
The hyJprid Monte Carlo will be described in next section. Here U is the energy of the cluster. The Rosenbluth 
factoiE^I (the normalization factor) for this step is calculated and denoted as wi. 

3. Generate a monomeric Si04 configuration from the Boltzmann distribution, under the constraint that two of 
the O atoms are feed at the positions of the pair of O'^' atoms. The Fixman potential associated with this 
constraint is zeroH3 Generation of the configuration is done by a short hybrid Monte Carlo run. We show in 
Appendix ^ that the probability for generating this configuration is 

/Q.^ X 47rroo^w(rsi,ro'')exp(-/3f7m) 

p(Si04,roo) = 7 rr^ , (10) 

p(roo)^i,4 



where U^a represents the energy of the Si04 monomer. 



4. Take the Si02 configuration from the monomer by dropping the two fixed O atoms, and graft this group onto 
the O^^' atoms in step 2 in a way so that the O atoms left in the monomer would have fallen on top of the 
O^^-* atoms. Multiple orientations of the fragment are generated by rotation of the Si02 group around the axis 
passing through the two O'-^-' atoms bonded to the Si atom. One of the orientations is chosen from a biased 
probability proportional to exp(— /3C/), U being the incremental energy associated with the Si02 group. The 
Rosenbluth factor for this step is calculated and denoted as W2- 

A Si02 removal move proceeds as follows: 

1. Pick a Si^^-' atom from the cluster at random. Remove this Si atom and the two bonded O*^^^ atoms from the 
cluster. 

2. Make local moves of the two remaining groups of atoms in the same way described in step ^ for the Si02 
insertion, except that p{r oo) is not included in the biased probability. Calculate the Rosenbluth factor wi. 

For both the Si02 insertion and removal move, a reverse move must be carried out so as to satisfy detailed balance.0 
The acceptance probability ratio for these moves can be derived from the partition function and the probabilities of 
proposing the moves. As shown in Appendix the acceptance probability is given by 

acc[0nsi,no n„s^ + i^„o+2] _ _2 »(00)^°^ ZiAW^^^/k2 Wi"'' 

acc[n„sj+i,„o+2 o„s,,„o] ^' ° „ (^Si^^^) '"^ exp(-/3i7m) " 

Here n(00) is the number of qualified O*-^^ pairs for the insertion move, n ^Si^^^^ is the number of qualified Si*-^-* 
atoms for the removal move, and fc2 is the number of trial moves for the Si02 insertion. 
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B. Hybrid Monte Carlo 

The hybrid Monte Carlo schemeEi has the basic idea that one can use molecular dynamics to generate global trial 
Monte Carlo moves. It has the advantage that a time step too big for molecular dynamics can be used without causing 
numerical instability. The molecular dynamicjs-|algorithm used to generate trial moves must be area preserving and 
time reversible to guarantee detailed balances Our hybrid Monte Carlo move equilibrates mainly the strong bond 
vibration and bending motions with short time scales. The time step for the molecular dynamics in hybrid Monte 
Carlo is chosen to be around 5 fs. To initialize each move, a Gaussian velocity distribution for the atoms is generated. 
The molecular dynamics then proceeds for a fixed number of time steps, and a Metropolis criterion is used to accept 
or reject such a move. 



C. CBMC Moves 



The configurational bias Monte Carlo method (CBMC), first developed by Frenkel, Smit, and de Pablo,t2lLa suc- 
cessfully samples complex energy landscapes by using local information when proposing moves. Our CBMC scheme 
proceeds by selecting randomly a Si-O*-^-* bond. If the cluster does not break into two independent clusters when 
this bond is removed, then we reject this move. This test is performed for successive bonds to be regrown. The two 
clusters so generated are reconnected by regrowth of the deleted fragments along the O'-'^-' — > Si direction unit by unit, 
beginning from the SiOs group bonded to the O^^). Each SiOg group is treated as a rigid unit in the regrowth. The 
reverse move is performed so as to satisfy detailed balance. This type of move equilibrates the Si-O-Si bending angles 
and the torsional angles around the Si-0 bonds. All the bond distances are preserved. 



D. Window Sampling 



Although the cluster distribution can be obtained in principle from a direct simulation in the grand canonical 
ensemble, the large free energy barrier to nucleation makes sampling of large clusters infrequent and statistically 
poor. The cluster distribution is best calculated by patching together the curves from several simulations sampling 
different ranges of nsi- This scheme is called window sampling .E3 A biasing potential of this form is added to the 
system so as to sample different clusters with roughly equal probability. We chose a linear function of nsi as a biasing 
potential, with the slope determined from a short trial Monte Carlo run. This bias was corrected for when computing 
the averages. 



V. RESULTS 



A. Modeling the Synthesis Condition 



We performed the cluster simulations in the grand canonical ensemble specified by zgi, zq, and T = 430 K, a typical 
temperature for zeolite synthesis. To estimate the fugacities of the Si and O atoms, we make use of the experimental 
data by the requirement that the equilibrium constant between monomers and dimers in simulation matches that of 
the synthesis condition. We assume that at pH — 12 the dominant monomeric and dimeric ions are Si(OII)202^ and 
812(011)403", respectively, and look for the equilibrium constant of this reaction 

2Si(OH)202" +H2O ^ Si2(OH)403" +2 0H-. (12) 



This equilibrium constant ispestimated from several reactions, ioduding the dimerization of neutral silicate species| 
the dissociation of silicates, and the dissociation of water .E3 The estimated value of the dimerization constant is 
10-3-2 M (mol/l). Since the dissociation constant of water at 430 K is about 10"" '^ M^, and the pH is 12, the OH" 
concentration is estimated to be iQ-n '^+i^ = 2.O M. Therefore, 

^ = ^ 10-- M-, (13) 

where Cd and Cm are the concentrations of the dimer and the monomer in units of M, respectively. On the other hand, 
from eq. || we can deduce 



Cd _ 1!2 4!2 Z2J 



2! 7! zoZ,/ 



(AVno.), (14) 
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where Cd and Cm are the concentrations of the dimcr and the monomer in units of (no. /A ), respectively. Using 
4 = exp(311.12) A"'^^ together with a short Monte Carlo run to determine the left hand side of eq. we find 

^2,7 — exp(623.00) A (see Appendix ^). Equating eq. and eq. |l^, we obtain zq — 0.465 A , which we use in 
all our simulations. Actually, simulations directly at this value of zq are rather difficult to equilibrate; we perform 
simulations at a somewhat higher zq and reweight the data to this value. The fugacity of Si is related to the reference 
system Si concentration by psi = zsi^O^-^i,4/(l'4!) (no. /A ). To estimate the Si concentration, we consider a typical 
silicate solution including 8% of Si02 and 92% of water. Assuming that about 30% of the silica dissolve in water and 
that the density of the solution is 1 (g/cm^), we find an estimated Si concentration of 0.4 M. The existence of a gel 
phase may cause a range of non-uniform local concentrations, and so we actually explore a range of Si concentrations. 

Window sampling of different cluster sizes, nsi: is performed from nsi = 1 to nsi = 200. Each window covers 4 or 5 
values of nsi- Each run starts with a configuration taken from the preceding run. Five thousand Monte Carlo cycles 
were performed for 1 < nsi < 50, and at least 10^ Monte Carlo cycles were performed for 50 < ngi < 200. Each cycle 
consists of 50 Monte Carlo moves. The probability distributions obtained from these windows were patched together 
to get the free energies of cluster formation. Simulations with different initial configurations were used as a check to 
ensure the system is well equilibrated. No hysteresis was observed within the range of this study. 

B. Structure Analysis 

ComroQn-experimentally identified silicate species such as 4-ring (nsi = 4) and double 4-ring (cubic octamer, 
7iSi = 8)Ejcj were found in our simulation at small cluster sizes. Figure ^ shows the average number of silicon atoms 
as a function of connectivities for different values of nsi- The 1-connected Si atoms are hardly observed for ngi > 25, 
indicating that the cluster is relatively compact. Figure ^ shows a snapshot of a nsi = 77 cluster, which is nearly 
spherical. 

The observed bond length of Si-0 in our simulation at 430 K is 1.67 ± 0.04 A. Figure |5| shows the probability of 
the nearest neighbor Si-Si distance for a large cluster with nsi — 194. Only the distances between 4-connected Si 
atoms were calculated. As a comparison, we computed the probabilities for the T-T distance, T being Si or AJ, for 
both the 28 all-Si zeolites and the 68 aluminosilicate zeolites in the International Zeolite Association database. El The 
simulation distribution is peaked around 3.05 A, while the distributions from known zeolite structures are peaked 
around 3.1 A. Interestingly, the distribution produced from our all-Si atomic model is closer to the distribution from 
all-Si zeolite structures. We also calculated the probability of the Si-Si-Si angle, shown in Fig. ^. Also shown in 
the same figure are the distributions for the 24 all-Si zeolite structures and the 68 aluminosilicate zeolite structures. 
The all-Si distribution displays a sharp peak around 110°, which is close to the value for a perfect tetrahedron. 
The distribution has a shoulder around 90°, which implies an abundance of 4- rings. The small peak near 60° in 
the simulation data comes from 3-rings, which are occasionally found in zeolite structures. The asymmetry is these 
distributions in due to steric repulsion for small ring sizes. 

We calculated the ring size distribution for ngi — 77 and 'T'Sir^ 200 by counting in the cluster the fundamental rings, 
i.e., those rings that cannot be divided into two smaller ones.EJ As shown in Fig. the ring size is mostly between 4 
and 8, consistent with the structures of Si02 compounds. Note that in crystalline zeolite structures there are several 
non-planar, fundamental rings that are not normally considered to be rings, e.g. there is a fundamental ring with 12 
silicons that encircles the sodalite cage in the zeolite SOD. In the histogram of crystalline zeolite rings, therefore, we 
have eliminated all fundamental rings that are not sufficiently planar 

Larger rings appear when the cluster grows, and for nsi = 200 we observed up to 11-rings. The smooth distribution 
of ring sizes suggests the cluster is amorphous. Interestingly, 3-rings appear only for the larger clusters. As shown in 
Fig. 1^, the 3-rings are mainly distributed near the surface of the cluster. This result is in accmd with experiments on 
amorphous silica, where 3-rings are observed only on the relatively flat surfaces of the silicaJ13 

C. Free Energy of Clusters 

Figure || shows the cluster free energy at a set of different Si monomer concentrations. Table || lists the estimated 
nucleation barrier and the corresponding critical cluster size. Clearly, the nucleation barrier and the critical cluster 
size increase as the Si monomer concentration decreases. For psi = 0.017 M the cluster free energy of cluster does 
not pass through a maximum within the range of this study. The positive slope suggests that the monomer phase is 
stable at this concentration. The cluster free energy becomes almost linear for ngi > 100 at this density, suggesting 
that the bulk energy term starts to dominate the cluster energy. 
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To investigate the relative stability of the clusters predicted by our forcefield, we calculate the Helmlioltz free energy 
at 430 K for both the a-quartz and the MFI framework structure using thermodynamic intjegrationEO The reference 
state is taken to be the Einstein crystal, whose free energy can be calculated analytically.cJ Our calculation shows 
crystalline quartz is more stable in solution than theJijIFI structure by 14 kcal/mol, while solution calorimetry results 
give an enthalpy difference of 1.5 kcal/mol at 298 °K.t3 This seems to indicate that our forcefield over-stabilizes quartz. 
Such a tendency has been observed in other forcefields for Si, O systemso Although quartz is the thermodynamically 
most stable structure in air at this temperature, this may not be true in the solution or for small crystallite sizes. 
Cations and structure-directing molgcules affect the surface energy and pore size distribution, and thus change the 
equilibrium distribution of clusters.QEJ Therefore, the composite formed by the structure-directing molecules and 
ZSM-5 may be more stable than quartz under aqueous synthesis conditions. From our simulation we did not observe 
quartz-like clusters for nsi < 200. To see if quartz clusters are less stable than the clusters observed from simulation, 
we generated a few small clusters from the crystal structure of quartz and calculated the free energy in solution by 
thermodynamic integration at psi = 0.33 M. There is no unique way of choosing the atoms that constitute the cluster, 
and we choose the atoms such that the cluster is nearly spherical and has no 1-connected Si atoms. As shown in 
Fig. H, the calculated free energy of quartz clusters are all above the curve corresponding to psi — 0.33 M. This 
suggests the formation of a quartz nuclei may require a higher activation energy at this condition. Scattering in the 
data comes mainly from our arbitrary choice of clusters. When the cluster is large enough so that the bulk term 
becomes dominant, we calculate that the free energy of quartz decreases with a slope of -12.8/per Si atom, which is 
calculated from the free energy of a quartz crystal and the fugacities of Si and O atoms. The estimated asymptotic 
slope of the amorphous clusters identified in the simulation at psi — 0.33 M is -2.8/per Si atom. Thus, we expect 
quartz to become relatively more stable for large clusters. 



VI. DISCUSSION 



Our analysis of the clusters suggests that more open, amorphous structures are favored at low values of ngj. The 
Ostwald ripening rule says that the condensed phase that forms during nucleation is notpHoeessarily the thermody- 
namically most stable phase, but rather the phase that has the lowest nucleation barrier .ta'Cj It is possible that the 
quartz phase has a much higher nucleation barrier and is therefore favored only at large values of ngj. 

It is expected that both the concentration and pH value have an influence on the critical cluster size and the 
nucleation barrier. As shown in Fig. ^ both the critical cluster size and the nucleation barrier decrease as the Si 
monomer concentration increases. The pH value affects the critical cluster size and the nucleation barrier through 
the oxygen chemical potential, fxo- This effect can be investigated by reweighting the cluster distribution at different 
values of po, as shown in Fig. ^ Here the Si monomer concentration is chosen to be 0.33 M in all distributions. The 
sensitivity of the nucleation barrier to the pH value is characteristic of the sensitivity of zeolite synthesis to solution 
conditions. A decrease of pH would shift the equilibrium of eq. |l^ to the right and lead to lower nucleation barriers 
and smaller critical clusters. This does not mean that lowering the pH value always leads to faster crystallization, 
since other factors such as insufficient supply of monomers may prevent the formation of the nucleus of the crystal. 
Indeed, the mechanism of nucleation for zeolites seems to operate only at high pH, where the critical cluster size is 
on the order of 50 silicons. 

We make a rough estimatejif the nucleation rate in this system. The steady state nucleation rate, P, is given from 
classical nucleation theory aa£3 



\ bmi* J 

Here D denotes the diffusion coefhcient of the monomers in the solution, A stands for the atomic jump distance 
in the solution, n* is the critical cluster size, Ap is the free energy difference between the monomer phase and 
the condensed phase, Afl* is the nucleation barrier, and Na is Avagadro's number. We pick psi = 0.33 M and 
use the simulation results n* — 32, /3|A/i| = 2.8, and /3Af2* = 91.6. The diffusion coefficient D is estimated to 
be 10~^ cm^/s, and the atomic jump distance A is estimated to be 1 A. Substituting these values into eq. |l5|, we 
obtain P = 5.5 x 10~^^ (s~^-A~^). We calculate the rate of nucleation events under typical non-templated synthesis 
conditions of a solution volume of about 1000 cm^ and roughly 1 hour. The probability of nucleation predicted for 
such an experiment is 1. 9.. This rough estimate is closer than is typical for comparisons of homogeneous nucleation 
theory and experiment .cZlc3 While silicate nucleation is often thought to be a strongly heterogeneous process, the 
present results suggest that the degree of heterogeneity may not be as strong as commonly thought. Indeed, many 
zeolite species can be synthesized without the use of a structure directing agent, and it may be that local fluctuations 
in solution conditions can be responsible for enhancing the rate of nucleation. Even under the assumption of a 
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perfectly homogeneous nucleation event, discrepancies between our estimate and the true rate may result from use of 
a simplified forcefield, uncertainties in the nucleation barrier, and assumptions in classical nucleation theory. 

Large rings with more than 10 constituent Si atoms are found in a few zeolites. These rings form one or two 
dimensional pores that are responsible for unique adsorption and catalytic selectivities. However, synthesis of zeolites 
with large rings is difficult and almost always requires the use of structure directing molecules such as tetrapropyl 
ammonium ions. The low population of large rings in the clusters observed in the present simulations, therefore, can 
be attributed partly to the absence of structure directing molecules. 

The synthesis condition modeled in this work strongly favors clusters with higher connectivities. It is observed 
during the simulation that Si and O atoms on the surface rearrange themselves to form bonds and thus minimize the 
cluster energy. As a result, 1-connected Si atoms are hardly observed for clusters larger than 25 Si atoms, as shown 
m Fig. |. For very large clusters we expect that almost all Si atoms are 4-connected. 

The nearest neighbor Si-Si distance observed in the cluster, when compared with the distribution calculated from 
zeolites, tends to be more condensed. Nevertheless, it displays a broad distribution, and therefore large distances that 
can lead to expanded structures are also found. The probability distribution for the Si-Si-Si angle observed in the 
cluster is closer to the distribution of the 68 aluminosilicate zeolites. Aluminum atoms are known to increase disorder 
in the crystal structures. For example, highly strained 3-rings are found only in zeolites having aluminum or other 
non-silicon elements. Such disorder is also found in the finite size clusters studied here, as indicated by the presence 
of smaller rings such as 3-rings and 4- rings. Note that the experimental peaks corresponding to 4-rings and perfect 
tetrahedral order are reproduced in the simulation. 



VII. CONCLUSION 



We have constructed an atomic-scale model for the silicate solution during zeolite synthesis in the absence of 
structure directing agents. Novel Monte Carlo moves were developed to sample the equilibrium distribution of silicate 
clusters. The nucleation barrier is estimated to be on the order of 10^ ksT, and the critical cluster size is estimated 
to be between 25 and 50 silicons. Interestingly, the nucleation behavior is rather sensitive to the Si concentration in 
and the pH of the solution, in accord with experimental experience. Structural analysis shows the results reproduce 
physical properties of condensed phases of silica, such as the distance distribution and topological features. The 
smooth distribution of ring sizes suggests the cluster is amorphous up to ngj — 200. Future work is planned to take 
into account the effect of ions and structure-directing molecules on the energetics of the clusters. 
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APPENDIX A: GENERATING A MONOMER CONFIGURATION WITH A FIXED roo 

We consider the configurational integral of a monomer with two of its O atoms fixed. Let roo be the vector that 
spans the two fixed O positions. Let p(roo) droo be the probability that two chosen O atoms of a monomer have a 
distance between roo and roo + droo- Then 

piroo)droo = ^"""^"f /drsidro^<5(Rcn.)«^(rsi"^', ro"°) 

xS{ro2 -roi - roo)exp[~/?[/(rsi,ro'*)], (Al) 

where the integral at the right hand side is the constrained configurational integral of our interest. Therefore, the 
probability for generating the constrained monomer configuration, as given in cq. is 

/c-n \ w(rsi,ro*)exp(-/3;7) 
p(Si04,roo) - 



J drsidro'^S(Rcm)S{ro2 - i"Oi - roo)w(rsi, ro'') exp[-/3t/(rsi, ro"*)] 
47rroo^w(rsi, ro"^) exp{-(3U) 
Piroo)Zi,4 
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The probability p(roo) is predetermined with a short hybrid Monte Carlo run. 

APPENDIX B: THE ACCEPTANCE RATIO IN GRAND CANONICAL MONTE CARLO 

Consider a simple, monatomic gas in a grand canonical ensemble. The grand canonical partition function of the 
system is 

^ - /dr^exp[-/3C/(r^)], (Bl) 

N ' '' 

where z = exp(/3^)/A^ is the fugacity, N is the number of particles, f3 is the reciprocal temperature, and U is the 
energy of the system. The Nl comes from the indistinguishability of the particles, and it is used here t o re move 



overcounting of the states. A more convenient way to think about the density of the states is to rewrite eq. Bl 



2 = ^ / d{v^} exp [-PU ({r^})] , (B2) 

where {r^} is the set of coordinates without individual labeling of atoms. Thus, the probability of a state is given by 

Piiv"}) = ^exp[-/3C/({r^})]. (B3) 

The detailed balance condition says we must accept the particle insertion and particle removal moves by the following 
criterion: 

acc({r^-}^{r^+i}) ^ pi{r^+^}) a{{r^+^} ^ {r^}) 

acc({r^+i} ^ {r^}) p({r^}) a({r^} ^ {r^+i}) ' ^ ' 

where a denotes probability of proposing a move, and acc denotes probability of accepting the proposed move. Let p 
be the probability of attempting a particle insertion or a particle removal move. We choose to perform an insertion 
or a removal move with equal probability. For the insertion move we insert a particle in the system volume V at 
random. Thus, 

a ({r^} - {r^+i}) ^ p^i (B5) 
For the removal move we choose at random a particle to be removed, and 

a({r^+i} - {r^}) = P^j^, (B6) 



Substituting eqs. B5 and |B6| into eq. B4, we get 



This is the usual acceptance ratio. E3c3 

In the present case, we have a two component system, so for the fixed center of mass single cluster 

p({rsi"-'},{ro"°}) = ° exp [-/?[/ ({rsi"^-},{ro"°})]. (B8) 

I — ^cluster 

For insertion or removal of a Si02 fragment, detailed balance reads 

acc[o„s.,no ^ n„s.+i,„o+2] ^ p({rsi"si+i},{ro"o+^}) 
acc[n„si+i,„o+2 o„si,„o] p ({rsi"si}, {ro"o}) 

a [n„gj+i^„o+2 o„sj,„o] 



a [Onsi,no ~^ n„sj + i,„o+2j 



(B9) 
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Here Onsi.no represents the old state with particle positions ({rgi"^'}, {ro"°}), and n„g;+i_„Q+2 represents the new 
state with particle positions ({rsi"^'"''^}, {tq"'-'^^}). Let ki be the number of trial moves generated for the two Si 
groups, and let k2 be the number of trial moves generated for the Si02 fragment. We choose ki = 50 and ^2 = 50. 
The probability of generating a Si02 insertion move, a[Onsi,no ~^ "'^nsi+i,no+2], is the product of probabilities of each 
selection process as described in the text. Trial moves for the Si02 removal move must be generated so as to calculate 
the acceptance probability. Detailed balance in this case describes the transition probabilities between two states 
that are not only characterized by their particle coordinates but also by the trial moves generated in the intermediate 
steps. We impose the super-detailed balance conditionj^j It follows that 



1 1 ki\ p{roo)exp(^-pUi^''' 



2!exp(-/3C/2(")) 
x27rp(Si04,roo) ^ ^ 

(fci - 1)! 

fel-1 



(27r)'=2 



1 w(rsi,ro'^) 
i(00)(°) Zi,4 



exp 



-/3(t/„ + C/i(")+C/2("))' 



(n) (n) 

fci!(fci-l)! fca! 
■ Ci^'^i-i (27r)fe2-i^ 



(BIO) 



JJi^") and /72^"^ denote the energy of the two chosen Si groups and the inserted Si02 group, respectively, and Ci is a 
constant that is related to the phase space volume generating trial moves for the two Si groups. The value of Ci is 
irrelevant because it cancels in the acceptance ratio. The factor of 1/2 accounts for the two possible choices for labeling 
groups 1 and 2. The factor 2tt before p(Si04,roo) accounts for the fact the rotational angle of the Si02 fragment 
has not been chosen yet; equivalently, this factor accounts for the values of the rotational angle that contribute to the 
flux. The factor of 2! accounts for the number of indistinguishable Si02 configurations contributing to the flux. The 

Rosenbluth factors are defined by wj"^ = X^tiLo P(''oo exp ^—/3f/;^'''^^ /(47rroo^) and W2 — '}2^k=o^'^'P(~P^'t^)- 
Note that a biased insertion volume instead of the system volume is used. The fci! and ^2! account for the number 
of indistinguis hable trial moves that contribute to the flux. The term (ki — 1)\/Ci'^~^ at the right end of the first 
equality of eq. BIO accounts for the probability of generating fci — 1 trial moves for the corresponding Si02 removal 
move. These fci — 1 trial moves combined with the configuration chosen in the insertion move are used to calculate 



the Rosenbluth factor w^°^ , which is described below. Inserting eq. A2 into the right hand yields the second equality. 
Since we are considering transitions between accessible configurations, the geometrical cutoff function w(rsi,ro^) is 
always one and can be omitted from eq. B10| . 

The probability of generating a Si02 removal move is 



1 1 fci! exp(-/3C/iM 



/ (n) n (o) 

(fci - 1)! (fc2 - 1)! 
^ Ci^^-^ (27r)'=2-i 



(BU) 



where w\"' = Yl'k=o^'^v{~l3U'^^). Note that wf^' and w^^' are asymmetric because in the Si02 removal move 



4°'' = EfcLoexp(-/3;7}''^). Note that wf^ and 
p(roo)/(47rroo^) is not used to bias the selection process. The first term at the right hand side accounts for the 
probability that we choose the appropriate 2-connected Si atom at random. We always choose the two O^^^ atoms 
bonded to the chosen Si'^' atom. The factor of 1/2 accounts for the two possible choices for labeling groups 1 and 2. 
Substituting eqs. B8, BIO, and Bll into eq. B9, and using C/2''"-' 



thus complete the proof. 



C/i'") - Ui^°^ = C/(") - U^°\ we obtain eq. |rT| and 
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APPENDIX C: THE REACTIVE MOVES 

In this section we describe the other Monte Carlo reactive moves. A SiOa insertion move proceeds as follows: 

1. Pick an O*^^^ atom from the cluster at random. 

2. Generate a monomeric Si04 configuration from the Boltzmann distribution at the desired temperature. This is 
done by a short hybrid Monte Carlo run. 

3. Take the SiOa configuration from the monomer by dropping the last O atom, and graft this group onto the 
O^^^ in step 1 in a way so that the O atom left behind would fall on top of the the O'^-*. The orientation of 
the fragment is chosen from a biased probability of the cosine of the new Si-O-Si angle and the two associated 
torsional angles. Within the cluster and monomer partition functions, the two associated torsional angles and 
the orientation of the monomer are already randomized. We optionally choose to re-randomize these. The 
biased probability is chosen to be exp(— /JVaj/C, where is given in eq. ^ for the O atom and 

C = j d{cose)exp{-pV3) (CI) 

A SiOa removal move picks a Si^^-' atom at random and removes the Si^^-' atom and the three associated O'-^-' atoms 
from the cluster. The acceptance probability is given below: 

acc[0ns,,„o nn3^ + i,„o+3] _ ^ 3 (O^^^)^"'' C(") Zi^4 

acc[n„3;+i^„o+3 ^ o„si,„o] ' n^^Si^^^y"^ ^ 3! exp(-/3C/m) 

xexp[-/3(AC/-y3^"^)], (C2) 

where AU is the change of system energy, and Um is the energy of the monomer. The factor of 1/2 accounts for the 
fact that the monomer orientation is already random in the absence of a bias; equivalently, this factor accounts for 
the values of cos 9 that contribute to the flux. The factor of 3! accounts for the number of indistinguishable Si03 
configurations that contribute to the flux. 

The idea of the hydrolysis/condensation move is to couple an O atom insertion/removal move to a local re-orientation 
of the nearest atoms, so as to alleviate part of the energy barrier present along the reaction path. We demonstrate in 
Fig. ^ how a typical hydrolysis/condensation move alters the local configuration of a cluster. The hydrolysis move 
proceeds as follows: 

1. Pick randomly a Si-O*-^-* bond that is to be broken later. If when this bond is removed, the cluster breaks into 
two independent clusters, reject this move. Otherwise, mark the two Si atoms to which this O*^^^ is bonded. The 
"first" Si atom is the one in the chosen Si-O'^^' bond. We name the O''^-' atom in the Si-O'^^-' bond the pivot O 
atom. Count the O^^-* connectivities of the two Si atoms after the Si-0 bond has been broken. Define the two 
connectivities to be ii and 12, respectively. Denote as group 1 the first Si and the (4 — ii) O'^^^ atoms bonded 
to it. Denote as group 2 the second Si and the (4 — 12) O*-^-' atoms bonded to it. 

2. We seek possible low energy positions for the O atom to be inserted by defining an insertion volume around the 
first Si atom. The way this volume is defined is quite arbitrary. This insertion bias affects the efficiency of the 
moves but does not alter the equilibrium results. We pick the insertion volume Vins according to ii: 

• If ii — 3, the insertion volume is defined as a sphere centered at the pivot O atom with a radius Re- 
Therefore, 

Vins = ^^i?c'- (C3) 

• If ii = 2, the insertion volume is a torus defined by revolving a circle around an axis passing through the 
two O*-^^ atoms bonded to the first Si atom. The circle resides on the plane defined by the two O*^^) atoms 
and the pivot O atom, is centered at the pivot O atom, and has a radius Rc- We find that for this choice 

Vins = 2^2^,evi?c', (C4) 

where Tj-ev is the distance between the pivot O atom and the axis of revolution. 
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• If zi = 1, the volume is defined as a shell formed by two concentric spheres centered at the only O*^^^ atom 
bonded to the first Si atom, each having a radius of + R^. and — Rc, respectively. Here is the 
distance between the O'^^-' atom and the pivot O atom. The volume is 

Mns = ^7ri?e(3rd2 + i?e2). (C5) 

Note that the radius i?c can be chosen independently in different cases. We use R^ = 0.3 A in all cases. Now, 
insert the new O atom randomly in the volume Vins, and move group 1 as a rigid body, so that the Si atom in 
group 1 moves away from the pivot O atom and bonds to the inserted O atom. This move is done according to 
the position of the insertion O atom: 

• If ii = 3, do not move group 1. 

• If ii = 2, rotate group 1 around an axis passing through the two O'-^-' atoms bonded to the Si atom, so 
that the plane defined by the two O*^^^ atoms and the pivot O atom, when rotated, will coincide with the 
plane defined by the two O*^^^ atoms and the inserted O atom. 

• If ii = 1, rotate group 1 so that the axis passing through the O*^^^ atom and the pivot O atom, when rotated, 
will coincide with the axis passing through the O*^^^ atom and the inserted O atom. This is followed by a 
random rotation of group 1 around the axis. 

3. Next, move group 2 according to the connectivity i2'. 

• If 12 = 3, do not move the Si atom. Choose randomly a new position for the pivot O atom from a sphere 
centered at the pivot O atom and having a radius Rc- 

• If 12 = 2, rotate group 2 by a random angle around the axis passing through the two O*-^-* atoms bonded 
to the second Si atom. 

• If 12 = 1, choose a new orientation for group 2 from a uniform distribution of the cosine of the Si-O^^-'-Si 
bond angle and the two neighboring torsional angles. 

4. Configurational bias Monte Carlo is used in step || and step § to enhance the acceptance probability. The trial 
positions for the inserted O and trial moves for group 2 are generated k times, where k is an positive integer. 
We use k = 100. Each set of positions is picked with a probability proportional tot-the exp(— /?[/), U being the 
energy of the configuration. Calculate the Rosenbluth factor for these trial moves.E^I 

A condensation move attempts to remove an O^^^ atom and alter the cluster configuration so that the Si atom that 
lost an O atom becomes bonded to another O atom. The condensation move proceeds as follows: 

1. Pick randomly a pair of O'^-* atoms from a condensation list that is constantly updated during the simulation. 
The two O^^-' atoms must not be bonded to the same Si atom, and the two Si atoms to which the two O^^-* 
atoms are bonded must not be bonded to the same O^^^ The latter restriction avoids cyclic 2-Si rings. The 
two O atoms must also satisfy certain geometrical criteria that are used to eliminate geometrically impossible 
pairs so as to improve efhciency of the move. However, the list must include all pairs that can be reached by the 
reverse move, i.e., the hydrolysis move. In the condensation list the O atom pairs and are counted 
as two different pairs, while in the Si02 insertion list the O atoms are counted as one pair. Denote the Si atom 
bonded to the first O atom and the associated O'^^-' atoms as group 1, and denote the Si atom bonded to the 
second O atom and the associated O*^^-' atoms as group 2. Name the second O'"'^' atom from the list the pivot 
atom. 

2. Move group 2 atoms in the same way as described in step |^ of the hydrolysis move. 

3. Remove the first O^^^ atom from group 1. Reorient group 1 so as to form a new bond between the Si atom 
of group 1 and the pivot atom. This reorientation depends on the connectivity and is the same as step || of 
the hydrolysis move, except that the roles of the pivot atom and the inserted atom are reversed. Multiple trial 
moves are generated, and the Rosenbluth factor for these trial moves are calculated. 



The acceptance probability is 



acc[0ns.,„o ^ n„3.,„o+i] ^ ^ n(SiO)<°^ (qq-. 
no] n(00)(") 
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where n,(SiO) and r7,(00) are the number of quahfied Si-0 and 0-0 pairs, respectively, and w^'^^ and w^""^ represent 
the Rosenbluth factors for the forward and reverse moves, respectively. 

A prune-and-graft move attempts to "prune" a fragment off the cluster and "graft" it onto a different site of the 
cluster. It proceeds as follows: 

1. Pick a Si-O'-^-' bond at random. If the cluster stays a cluster when this bond is removed, then we reject this 
move because the branch so generated is not free. Otherwise, mark the entire fragment that is still connected 
to this O^^^ after cleavage of the Si-0^^^ bond. 

2. Pick randomly an O'^' that is not part of the marked fragment. 

3. Move the marked fragment in a way for the Si atom liberated in step 1 to bond to the O*^^^ chosen in step 2. 
During the move, the fragment is treated as a rigid body. Require that the bond has the same length as the 
bond broken in step 1. The orientation of the fragment is chosen from a biased distribution of the new Si-O-Si 
angle and from a uniform distribution of the two neighboring torsional angles. The biasing probability is the 



(C7) 



same as that used in SiOa insertion move, with the normalization constant C given by eq. CI 
The acceptance probability is 

acc(o,igj^„Q — s- n„3j^„Q 

acc(n„sj^„o o„si,„o) C(°) exp(-/3t/(°)) ' 

APPENDIX D: COMPUTING THE CONFIGURATIONAL INTEGRAL OF A MONOMER 

The configurational integral for a monomer, Zi^i, is calculated by thermodynamic integration from infinite temper- 
ature to the synthesis temperature: 

lnZi,4 = lnZi,4(/3 = 0)- / {U^)p'dl3'. (Dl) 

Jo 

Here {U„i)p' is the canonical ensemble average of the monomer energy at /?', and 

Zi4(/3 = 0) = y"drsidro(J(Rcm)w(rsi,ro'') 

= j (ixidx2rfX3C?X4(iRcmW(0, X^)(5(Rcm) 

4 „ 



3^i?max'j , (D2) 

where x^ = roi — rsii a-nd Rinayi = 2.6 A is the maximum Si-0 bond length. We compute the value of this integral by 
splitting the range of integration int o three regions, and we use 16-point Gauss-Legendre integration in each region. 



Performing the integration of eq. Dl, we find Zi,4 = exp(311.31) A^^ at T = 430 K. An alternative way of computing 
the free energy is to couple the cluster to a set of non-interacting harmonic oscillators whose free energy is known, using 
a parameter A to adjust the relative contributions of the two systems. Integration from A = to A = 1 gives the free 
energy difference between the monomer and the harmonic oscillators and, thus, the absolute free energy of the cluster. 
To compute the integral, we again split the range of integration into three regions and use 16-point Gauss-Legendre 
integration. The value obtained in this way agrees with the above result within statistical error, lnZi^4 = 311.12. 
Note that in the harmonic oscillator reference system, the atoms are distinguishable, and so an "identity exchange" 
move must be added to the simulation so that the free energy difference between the indistinguishable cluster and 
the distinguishable reference system is properly calculated. This distinguishability factor, as well as the rotational 
entropy, shows up as a large contribution at small A in the integral for the free energy difference. We adopt the 
Einstein crystal result, lnZi.4 = 311.12, as it is likely more acc ura te because the thermodynamic integration that 



gives this result avoids the singularity that exists at /3 = in eq. ^ . 

To test detailed balance for our SiOs insertion and removal moves, we compared the values of ln(Z2, 7/^1,4) — 311.88 
obtained from our simulation using only the SiOs insertion and removal moves to change the number of atoms with 
that from the integration via harmonic oscillators, ln(Z2,7/.^i,4) — 311.77. The results agree within statistical errors. 
To test detailed balance for the hydrolosis and condensation moves, we compared the free energy differences between 
a liner and cyclic tri-silicate species. We find the result ln(Z3_io/2^3.9) = 6.32 from a simulation employing only 
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the hydrolosis and condensation moves to change the number of atoms, and we find ln(Z3,io/^3,9) = 6.38 from the 
integration via harmonic oscillators. The level of agreement is, again, within statistical errors. To check the Si02 
insertion and removal move, we consider the equilibration between a dimer and a cyclic trimer. From the simulation, we 
find ln{Z3^g/Z2,7) = 307.84. From the integration via harmonic oscillators we find \n.{Z3^g/Z2,7) = 307.65. Again, the 
results agree within statistical errors. Finally, we check a combination of the hydrolosis and condensation moves and 
the Si02 insertion and removal moves. From the above thermodynamic integration results, we find ln(Z3_io/.^2,7) = 
314.03. From simulation, we find ln(.^3,io/.^2,7) = 313.96, which is again within statistical errors. 
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TABLE I: Constants in the interaction potential for silicate ions. The units of length and energy are A and e^lk= 14.39 eV, 
respectively. 





q(e) 


a (A3) 






Si 


1.60 


0.00 






o 


-0.80 (for O*^^) 


2.40 








-0.90 (for O'^^) 










•n 


H 






Si-Si 


11 


0.057 






Si-0 


9 


11.387 






0-0 


7 


51.692 








B 


/ 


(deg) 




Si-O-Si 


1.40 


1.0 


141.00 


2.60 


O-Si-O 


0.35 


1.0 


109.47 


2.60 



TABLE II: The nucleation barrier and critical cluster size at different Si concentration and pH — 12. At 0.017 M, the free 

energy of cluster does not pass a maximum up to nsi = 200; thus the nucleation barrier and critical cluster are undefined. 

Si coucoutratiou (M) ,JA(2 ((liiiicusiinik'ss) C'ritical Cluster (resi) 
0.75 66.2 27 

0.33 91.6 32 

0.08 126.8 45 

0.017 
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FIG. 1: The reactive moves change the number of atoms, as shown on the nsi, no lattice. The two dashed lines, no = 3nsi + 1 
and no = 2nsi, represent the upper and lower bounds for the possible lattice points. 



FIG. 2: A Si02 insertion/removal move. The upward arrow represents an insertion move, while the downward arrow represents 
a removal move. 



FIG. 3: The connectivities of Si atoms for different values of nsi. The statistical errors are smaller than the size of the points. 



FIG. 4: A snapshot of a cluster with nsi = 77. 



FIG. 5: The probability distributions for the adjacent T-T atom distances, T being Si or Al. Note that the cluster is simulated 
at T = 430 K, while the other two distributions are computed from crystal structures. 



FIG. 6: The probability distributions for the T-T-T angles, T being Si or Al. 



FIG. 7: a) Average ring size distribution at nsi = 77 and ngi = 200. b) Ring size distribution of known zeolite structures. 



FIG. 8: Probability of observing the ccntroid of a 3-ring at a given radius from the center of the cluster for a 197 Si-atom 
cluster. The distance is scaled by the radius of gyration of the cluster, Rg. 



FIG. 9: The free energy of cluster formation along the nucleation coordinate at pH = 12. The unit of the Si monomer 
concentration, psi, is M (mol/Z). The dashed line represents the free energy of cluster formation at psi = 0.33 M, beyond the 
critical cluster size. The free energies of quartz cluster formation are calculated using values of /usi and //o corresponding to 
psi = 0.33 M and pH = 12. The statistical errors are roughly the size of the points. 



FIG. 10: The free energy of cluster formation along the nucleation coordinate at different values of the pH and a fixed silicon 
monomer concentration of psi = 0.33 M. 



FIG. 11: A typical hydrolysis/condensation reactive move. The downward arrow represents a hydrolysis move, while the upward 
arrow represents a condensation move. 
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